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. ABSTRACT 

Context. The WMAP satellite has made available high quality maps of the sky in five frequency bands ranging from 22 to 94 GHz, 
with the main scientific objective of studying the anisotropies of the Cosmic Microwave Background (CMB). These maps, however, 
contain a mixture of emissions from various astrophysical origins, superimposed on CMB emission. 

Aims. The objective of the present work is to make a high resolution CMB map in which contamination by such galactic and extra- 
galactic foregrounds, as well as by instrumental noise, is as low as possible. 

Methods. The method used is an implementation of a constrained linear combination of the channels with minimum error variance, 
and of Wiener filtering, on a frame of spherical wavelets called needlets, allowing localised filtering in both pixel space and harmonic 
space. 

Results. We obtain a low contamination low noise CMB map at the resolution of the WMAP W channel, which can be used for a 
range of scientific studies. We obtain also a Wiener-filtered version with minimal integrated error. 

Conclusions. The resulting CMB maps offer significantly better rejection of galactic foregrounds than previous CMB maps from 
WMAP data. They can be considered as the most precise full-sky CMB temperature maps to-date. 
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t"~~ 1. Introduction fraction of the sky. For this reason, many ground-based and 

balloon-borne experiments have concentrated their observations 

OO The WMAP satellite is one of the most successful expen- in „ dean „ ions of ^ sk where ^fe. emission is low 
ments dedicated to mapping the Cosmic Microwave Background h tQ { neg ligibly the observations. For power spec- 
> (CMB). The all-sky maps obtained in the WMAP five frequency tmm estimation from full _ sk observa tions, a safe approach con- 
bands, in temperature and polarisation, offer the best data set to- sists in maski ions at low ^fc. latitude, and estimating 
date for making a sensitive all-sky map of the CMB anisotropies. pQwer spectra Qn ^ deanest regions of ±e sky The impact of 
H The CMB, however, is not the only source of emission at extragalactic point sources (evenly spread on the sky) on power 
WMAP frequencies. Diffuse galactic emission from several pro- spectmm estimates can be evaluated and corrected for using an- 
cesses contaminates the maps with an amplitude roughly propor- cm data (catalogues of known point sources and priors on the 
tional to the cosecant of the galactic latitude, compromising the statistical distribution of sources) 
observation of the CMB close to the galactic plane. In addition, 

a background of radio and infrared compact sources, galactic Besides the power spectrum, the CMB map itself is interest- 

or extra-galactic, contributes to the total emission even at high i"g for several additional purposes: 
galactic latitudes. Component separation consists in separating 
one or more of these sources of emissions from the others in the 

data. - As a CMB template, to be subtracted from millimetre-wave 
One of the main objectives of CMB experiments is the mea- observations when the scientific focus is on other emissions, 
surement of the CMB angular power spectrum Q which, with or to be used for the calibration of other instruments; 
the assumption of statistical isotropy, can be estimated on a - To assess the statistical isotropy of the CMB and check the 
homogeneity and isotropy of the Universe on the largest 
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To search for signatures o f non-trivial topology, as that of a 
multi-connected universe (Aurich et al. 2006; Caillerie et al. 
[2007)|NiSchou & Jaffe|20071 >; 
To search for correlations of the CMB map with other emis- 



sions (N olta et al.|2004{ |Fosalba & Ga ztanaga 2004 



Cabre 



et al.|2006llCao et al.|200 6, Pietro bon et al.|20061|McEwen 
et al.|2007[|Rassat et al.|2007| l; 

To search for signatures of non Gaussianity in the CMB 
(Mukherjee & Wang 20 041 |Vielva et al |2004[ |Cayon et al 



2005 , Cabella et al. 2006 ; McEwen et al. 2006b , 



Tojeiroetal. 



20061 |McEwen et a l. 2006a; Wia ux et al.||20071 |Yadav~& 



Wandelt|2007T i 



Several CMB maps obtained with the WMAP data are avail- 
able for such research projects. The WMAP team has released 
part-sky foreground-reduced maps in the Q, V and W bands, 
and maps obtained by an Internal Linear Combination (ILC) of 
all WMAP channels ( |Hinshaw et al.|2007[ l. |Tegmark et al.| ( |200"3"] l 
have produced CMB maps with WMAP one year data, and sub- 
sequently with three year data, based on an harmonic space ILC 
method. Eriks en et al.| ((2004 ) have an alternate version of the 
ILC CMB map at 1 degree resolution on one year data. Eriksen 
et al. (2007 1 use a Gibbs sampler to draw 100 realisations of the 
CMB under the posterior of a model of CMB and foregrounds; 
their estimated CMB map is the average of these realisations for 
three year data. On three year data again, [Park et al. ( 2007 1 use 
an ILC technique on 400 different pixel ensembles, selected by 
the spectral index of the foreground emission as estimated by the 
WMAP team using a Maximum Entropy Method (MEM). More 
authors have addressed component separation on WMAP three 
year data and produced versions of a "clean" CMB map (Maino 



et al. 2007 



Saha et al. 2007 Bonaldi et al. 2007 ). More recently. 



Kim et al.|(2008l have obtained a CMB map from WMAP five 



year data. 

All available maps suffer from limitations, some of which re- 
sult from specific choices in the way the CMB map is produced. 
Several of these maps, for instance, do not fully exploit the res- 
olution of the original observations. Some focus on cleaning the 
CMB from foregrounds at high galactic latitude, and are signifi- 
cantly contaminated by foregrounds in the galactic plane. Some 
are not full sky CMB maps. Finally, not all of the available maps 
have well characterised noise and effective beam. All these lim- 
itations impact their usefulness for accurate CMB science. 

In this paper, we address the problem of making a CMB tem- 
perature map which has the following properties : 



- being full sky; 

- being as close as possible to the true CMB (minimum vari- 
ance of the error) everywhere on the sky, and on all scales; 

- having the best resolution possible; 

- having well-characterised beam and noise. 



In the following, we first review these requirements and their 
impact on a CMB cleaning strategy (Section [2]). We then re- 
view and compare available maps in Section [3] In Section HJ we 
describe and explain our ILC needlet method. The approach is 
tested and validated on realistic simulated data sets (Section [5} 
before applying it to WMAP data (Section|6]l. We then compare 
our CMB maps to the other existing maps, discuss the results, 
and conclude. 

This paper considers only temperature maps. 



2. General considerations 

2.1. Requirements 

We start with a review of the requirements above, and on the 
implication on the method to be used. 

2.1.1. Full sky? 

WMAP data are full sky. We wish to devise a method which 
allows recovering an estimate of CMB emission everywhere, in- 
cluding in the galactic plane, and even (as much as possible) in 
the galactic centre as well as in pixels strongly contaminated by 
compact sources. 

The large scale correlation properties of the CMB make it 
possible to estimate the CMB emission even in unobserved re- 
gions, by some kind of interpolation. This is incidentally what 
is obtained with the Gibbs sampling technique of Eriksen et al. 
(2007). Equivalent in spirit although quite different in implemen- 
tation is the use of an "inpainting" method as that of Elad et al 



(2005). Such interpolation methods alone are not fully satisfac 
tory, as they discard information. In particular they do not allow 
recovering small scale CMB features in the mask. This is obvi- 
ous, for instance, in the Gibbs -sampling average map of Eriksen 
|eTaT1 ( |2007l l. 

At the opposite, one may try to separate components in the 
galactic plane independently of what is done at higher galac- 
tic latitudes, since levels and properties of foreground emission 
depend strongly on sky direction. |Hinshaw et al.| ((2007]) and 



Tegmark et al. (2003) divide the sky into several independent 
regions, perform component separation independently in these 
regions, and then make a composite map by stitching together 
these independent solutions. Such approaches discard informa- 
tion (zone-to-zone correlations) and require careful treatment at 
the zone borders to avoid discontinuities and ringing. 

A good method should perform well on both counts: lo- 
calised processing and full exploitation of large scale corre- 
lations of the CMB and of galactic foregrounds. This can be 
achieved with a spherical wavelet or needlet analysis of the 
maps (using the tools developed in Marinucci et al. (2008) and 
Guilloux et al. ( 2008) 1), which is our approach in the present 
work. 



2.1.2. Minimum variance? 

Recovery of a CMB map can be conducted following various ob- 
jectives, quantified by different "figures of merit". In this work, 
we choose, as most authors, to minimise the variance of the dif- 
ference between true and recovered CMB (this is the overall er- 
ror; it includes additive noise, foregrounds, and even multiplica- 
tive errors affecting the CMB itself). 

This choice alone does not fully characterise the method to 
be used. The best theoretical solution also depends on the model 
of the data. An overview of existing methods can be found in the 
review by Delabrouille & Cardoso (2007). 

In this paper, contrarily to other approaches which rely 
heavily on the structure of the data as described by a model, 
for instance a noisy linear mixture of independent components 
(Delabrouille et al. 2003), we assume as little as possible about 
the foregrounds and the noise. In fact, we assume nothing except 
the following: 



- The WMAP data are well calibrated with known beam in 
each channel; 
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The instrumental noise in all WMAP maps is close to be- 
ing Gaussian and uncorrelated; its pixel-dependent level is 
approximately known; 

The CMB anisotropies emission law is known to be the 
derivative with respect to temperature of a T=2.725K black- 
body; 

To first order, the template of CMB anisotropies is well rep- 
resented by a Gaussian stationary random field, the spectrum 
of which is given by the WMAP best fit (as will be seen later- 
on, this last assumption is needed only to derive the Wiener 
filter; it is not necessary for our needlet ILC map). 



2.1.5. Noise 

Throughout this paper, the term "noise" typically stands for 
all sources of additive error, i.e. instrumental noise and fore- 
grounds. 

2.2. Evaluation and comparison of reconstructed CMB 
temperature 

We briefly discuss here the tools used for characterising and 
comparing CMB temperature maps. 



These assumptions lead us to consider an "Internal Linear 
Combination" (ILC) method, followed by a Wiener filter to min- 
imise the error integrated over all scales. 



2.1.3. Best resolution? 

The WMAP data comes in five frequency channels with vary- 
ing resolution. To make the best of the data, we need a method 
which uses the smallest scale information from the W band, and 
information from additional bands (V, then Q, then Ka, and fi- 
nally K) at increasingly larger scales. Multi-scale tools are well 
suited for this purpose. 

The "best possible resolution" is not a well defined concept 
(and not necessarily the resolution of the W band). Indeed, there 
is a conflict between best possible resolution and minimum vari- 
ance, as one can smooth or deconvolve a map arbitrarily in har- 
monic space, reducing or increasing the total noise variance in 
the process. Here, we make a map at the resolution of the W 
channel over the full sky, leaving open the option to filter this 
map if needed to reduce the noise - or deconvolve it for better 
angular resolution. Note that additional global filtering or decon- 
volution does not change the signal to noise per mode (only the 
integrated S/N). 

The minimum variance map is obtained with a Wiener filter, 
which smoothes the map depending on the signal to noise ratio. 
As the noise and the contaminants are inhomogeneous, there is 
a strong motivation for the smoothing to depend on the location 
on the sky (the optimal solution to the resolution-variance trade- 
off depending on the contamination level, which is local). If we 
relax the constraint about beam homogeneity, again, spherical 
needlets offer a natural way to obtain such location-dependent 
smoothing. In the present, however, in order to preserve the con- 
stancy of the effective beam over the sky, we implement the 
Wiener filter in harmonic space. 



2.1.4. Accurate characterisation? 

A fully accurate characterisation of the beam and noise is not 
straightforward, in particular because of the limited knowledge 
about the original frequency maps, which automatically propa- 
gates into the final CMB map. This work makes several approx- 
imations about beams and noise. Beams are assumed symmetric 
and therefore described by the be transfer functions provided by 
the WMAP team. The instrumental noise is assumed uncorre- 
lated, although non stationary. Analytical analyses and Monte- 
Carlo simulations are used to characterise the residual noise of 
the final map, as well as to estimate the contribution of resid- 
ual foregrounds, and biases if any. This will be detailed further 
later-on. 



2.2.1. Map description 

A pixelised map is fully characterised by the specification of: 

- A set of temperature values y p in a number of pixels (here 
indexed by p); 

- The effective beam at each pixel p, which in the most general 



case is a function b 



p,p'> 



- The noise n p , the statistical properties of which, in the 
Gaussian case, are fully described by a noise covariance ma- 
trix Npj. 

The map value y p is then linked to the true signal value s p by: 

y P = 2_jbp, p >s p > +n p (1) 

p' 

The full characterisation of a given CMB map requires the 
specification of the additive noise n p and of the response b p ^ . 
When the beam is stationary over the sky and symmetric, which 
we assume in this work, it is fully specified by the coefficients be 
of the expansion of the beam in Legendre polynomials. 

2.2.2. Assumptions 

Throughout this paper, the beam is assumed symmetric. 
Although this is an approximation, most pixels of the WMAP 
map are "visited" by any particular detector through a wide 
range of intersecting scans. The average beam in that pixel then 
is an average of the physical beam over many orientations, which 
makes the symmetry assumption reasonable. 

In addition, in absence of any specific localised processing, 
the beam is assumed to be invariant over the sky. 

With the above two assumptions, the effect of beam convo- 
lution is best represented in harmonic space, with a multiplica- 
tive coefficient be, independent of m, applied to the harmonic 
coefficients cie m of the map. We assume perfect beam knowledge 
as well as perfect calibration, so that no multiplicative uncer- 
tainty is attached to the map description (the beam integral, ap- 
proximated as Yjp> b P , P ', is equal to unity independently of p, or, 
equivalently, the value of be f or I — is assumed to be exactly 
unity). 

The noise n p of the original WMAP maps, for each fre- 
quency channel and each differencing assembly, is assumed un- 
correlated from pixel to pixel, i.e. N PiP > - (n p n p >) = 5 pp >cr p . The 



variance cr is pixel dependent because of uneven sky coverage. 
Noise is non-stationary, but assumed to be Gaussian distributed. 



2.2.3. Comparison of maps at different resolution 

The comparison of CMB maps is meaningful only if the maps 
are at the same resolution. As long as beam transfer function 
does not vanish at any useful I (which is always the case for 
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Gaussian beams, and is also the case for WMAP best fit beams 
for all meaningful ranges of (), the resolution of any map can be 
changed to anything else by harmonic transform and multiplica- 
tion by the ratio of the beam transfer functions. This property is 
widely used throughout this paper. 

2.2.4. Masking 

We define "tapered" regions of the sky for map comparison at 
varying galactic latitude. In particular, we define a Low Galactic 
Latitude (LGL) region and a complementary High Galactic 
Latitude (HGL) region. The LGL region, used to evaluate results 
in the galactic plane, cuts completely all data above 30 degrees 
galactic latitude (and below -30 degrees), and has a 15 degree 
transition zone with a cosine square shape. All pixels at absolute 
galactic latitudes below 15 degrees are kept with a coefficient of 
1. The HGL region is the complementary, i.e. HGL = 1-LGL. 
These "tapered" regions allow the computation of local power 
spectra with negligible spectral leakage of large scale power into 
small scales. 



2.2.5. Power spectra comparison 

For a given beam (i.e. multiplicative response as a function of €), 
the comparison of the total map power as a function of € (i.e. of 
the power spectra of the maps), is a direct figure of merit. The 
lower the power spectrum, the better the map. 

Power spectra are computed independently for different re- 
gions of the sky (e.g. inside or outside the galactic plane). To 
minimise aliasing due to sharp cuts, we use masks with smooth 
transitions, defining the LGL and HGL regions described above. 

The power spectrum of a map is evaluated as: 



Ci 



1 



(2£+\)a 



2_j \ a tr>- 
m=-t 



(2) 



a is a normalising factor computed as the average value of the 
squared masking coefficientsPl 

Power spectra estimated directly in this way from a masked 
sky are unreliable for modes corresponding to angular sizes 
larger than the typical size of the zone of the sky retained for 
computation. 

2.3. Methods 

We now give a brief introduction to the two main methods used 
in this paper (ILC and Wiener filtering). Many other methods 
exist for CMB cleaning (or component separation in general), 
which assume varying degrees of prior knowledge about sky 
emission, and model the data in different ways. These meth- 
ods are not discussed nor used in this paper. For a review, see 
IDelabrouille & Cardoso|p007] l. 

2.3.1. The ILC 

The data are modelled as 



x — as + n 



(3) 



where x is the vector of observations (e.g. five maps), a is the 
response to the CMB for all observations (e.g. a vector with 5 



1 The masking coefficient is simply 1 in regions kept for power spec- 
trum computations, in regions masked, and between and 1 in the 
transitions. 



entries equal to 1 if WMAP data only are considered) and n is 
the noise. Here it is assumed that all observations are at the same 
resolution. 

The ILC provides an estimator Silc of s as follows: 



a'R- 



SlLC - 



(4) 



where R is the empirical covariance matrix of the observations 
(e.g., a 5 x 5 matrix when 5 channels are considered). 

Note that the ILC solution of Equation |4] is the linear filter 
which minimises the total variance of the output map, under the 
condition that the filter has unit response to the signal of interest 
(the signal with emission law given by vector a). 

The details of the method for its implementation in the con- 
text of this work are further discussed in Section [4] 

2.3.2. Wiener filtering 

Given a single CMB map of known beam (assumed to be con- 
stant over the sky), it is possible to minimise the contamination 
by noise and foregrounds by (one-dimensional) Wiener filtering. 
The data is modelled as x - s + n, where now x is a single map, 
s the true CMB and n the noise. The Wiener filter gives to in- 
dividual "modes" a weight proportional to the fraction of signal 
power in that mode, i.e. 



Hc e 



Stir. 



b 2 ,C l + Nt 



Xlil 



(5) 



where b 2 ( Q and N( are the power spectra of the (smoothed) CMB 
and of the noise (including smoothed foregrounds or foreground 
residuals) respectively, and xtm is the original noisy CMB map. 
It should be noted that if the CMB and the noise are uncorre- 
lated, then b 2 ( C( + Nt — Xt is the power spectrum of the map 
xtm, and the Wiener filter can be estimated directly using only a 
prior on the CMB power spectrum (assuming Q is known), and 
estimating Xf on the map itself. 

Wiener filtering in harmonic space minimises the variance of 
the error in the map if signal and noise are Gaussian and station- 
ary. For non-stationary contaminants, the Wiener filter (J5l) is still 
meaningful, but is no longer optimal. Some efficiency may be 
regained by an implementation in another domain than the har- 
monic space (e.g. needlets), but this is not investigated further in 
this paper. 

3. Evaluation and comparison of available maps 

Before describing how to make yet another CMB map from 
WMAP data, we review the existing maps and evaluate in which 
respect they can be improved. 

We start with a discussion of the existing methods, identify- 
ing for each of them strengths and weaknesses of the approach, 
and their foreseeable consequences. Available CMB maps ob- 
tained from WMAP data are compared in terms of resolution, of 
the estimated contamination by foregrounds and of noise level. 
In absence of an absolute reference, discrepancies between avail- 
able maps are also evaluated. This comparison permits to esti- 
mate typical uncertainties, and to outline the "difficult regions" 
for CMB cleaning (which, unsurprisingly, are mostly located 
close to the galactic plane). We also look specifically for residu- 
als of galactic contamination by comparing the power spectrum 
of the reconstructed CMB map at high and at low galactic lati- 
tudes. Significant discrepancies between the two are interpreted 
as indicative of a residue of foreground emission. 
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Table 1. Available CMB maps. 



NAME resolution data used Reference 

WILC1 P Pyr Bennett et al. } 2003 F 

TILC1 W channel 1-yr [legmark et aL](|:_UU3} 

EILC1 1° 1-yr briksen et al. 1 2004 1 



WILC3 


1° 


3-yr 


EGS3 


3° 


3-yr 


PILC3 


1° 


3-yr 


TILC3 


W channel 


3-yr 



Hinshawetal. (2007) 
iinksen et al. ( 2007 ) 
|rark et al.||Z007 } 



URL 



http://lambda.gsfc.nasa.gov/product/map/drl/imaps_ILC.cfm 

http://space.mit.edu/home/tegmark/wmap/cleaned_map.fits 

http://www.astro.uio.no/~hke/cmbdata/WMAF_lUJ_lagrange.hts 



1° 5-yr Gold etal.l 2008) 

1° 5-yr Kim et al. (2008) 



http://lambda.gsfc.nasa.gov/product/map/dr2/ilc_map_get.cfm 

http://www.astro.uio.no/~hke/gibbs_data/cmb_mean_stddev_WMAF3_nb4_3deg.hts 

http://newton.kias.re.kr/~parkc/CMB/SILC400/SILC400_bc.fits 

http://space.mit.edu/home/tegmark/wmap/cleaned3yr_map.fits 



WILC5 
KILC 



http://lambda.gsfc.nasa.gov/product/map/dr3/ilc_map_get.cfm 

http://www.nbi.dk/~jkim/hilc/ 



3.1. Available maps 

3.1.1. The WMAP ILC 

The ILC maps obtained by the WMAP team (denoted as WILC1, 
WILC3 and WILC5 hereafter, depending on whether they are 
obtained using one year, thre e year or five year da ta) a re de- 
scribed in|Bennett et al.| ( |2003] ), pir__"haw et al.| ( [2007] >, and [Go!dl 
et al. ( 2008| ) respectively. 

For the three year and five year maps, the original frequency 
maps are smoothed to a common resolution of one degree. The 
sky is subdivided into 12 regions. One big region covers most of 
the sky at moderate to high galactic latitudes. The rest of the sky, 
concentrated around the galactic plane, is divided into 1 1 regions 
of varying galactic emission properties (amplitude and colour). 
An Internal Linear Combination is performed independently in 
each of these regions. A full sky composite map is obtained by 
co-adding the maps of the individual regions (with a ^ 1 degree 
transition between the zones to avoid sharp edge effects). Finally, 
a bias (unavoidable consequence of empirical CMB -foreground 
correlation) is estimated by Monte-Carlo simulations, and sub- 
tracted from the composite map, to yield the final CMB map. 

The three year and five year ILC maps differ from the one 
year map in several respects. The most significant is the recog- 
nition of the existence of a bias, and the attempt at correcting 
it using simulations. The limitations of the maps include their 
resolution (one degree), and the use of small regions in the ILC, 
which is bound to cause more bias than necessary on large scales 
(comparable to patch sizes), as well as edge effects. This results 
in discontinuities between regions, obvious for instance in the 
estimated bias map shown in Hinshaw et al. (2007 1. 

It is worth mentioning that in the method used by the WMAP 
team, the coefficients of the linear combination used over most 
of the sky (region 0, which corresponds to the largest part of 
the sky at high galactic latitudes and a few low galactic latitude 
patches, and region 1, in the galactic plane but away from the 
galactic center) are set using only a small subset of the data in- 
side the Kp2 cut (where the galactic emission is the strongest). 
This choice favours the rejection of galactic contamination, at 
the price of sub-optimal weighting of the observations in regions 
where the error is dominated by noise. It also assumes that the 
emission laws and relative power of the different foregrounds 
are the same in these regions, which is a strong (and probably 
wrong) assumption. 

Furthermore, the ILC weights are set by minimising the vari- 
ance of the map at one degree resolution. Modes at higher I get 
very sub-optimal weighting, as they do not contribute signifi- 
cantly to the total variance of the one degree map. The K and 
Ka bands, in particular, contribute respectively to about 0.156 
and -0.086 for region (the largest one) for WILC3. As a con- 
sequence, the final ILC map can not be meaningfully decon- 



volved to better resolution than about 1 degree (because this 
would blow-up small scale noise coming from the lowest fre- 
quency channels). 

Finally, there is also an unsatisfactory degree of arbitrariness 
in the choice of the regions, which depend on priors about fore- 
ground emission, and are somewhat elongated across the galactic 
plane for no particular reason. Although none of these choices is 
unreasonable, the impact of this arbitrariness on the final map is 
difficult to evaluate. 

For all these reasons, the WILC maps leaves considerable 
margin for improvement. We aim, in particular, at obtaining a 
CMB map with better angular resolution, and with a better han- 
dling of non-stationarity and scale dependence of the contami- 
nation (foregrounds and noise). 

3.1.2. The WMAP foreground-reduced maps 

For temperature power spectrum analysis, the WMAP team 
has used part-sky foreground-reduced maps. The processing for 
foreground removal for the three year and five year releases is 
described in Hinshaw et al. (2007). Model templates for galactic 
emission are fitted to the Q, V and W WMAP channels outside 
of the Kp2 mask. A linear combination of synchrotron, free-free 
and dust, based on this fit, is then subtracted from the full sky Q, 
V, W maps. 

In this procedure, a first galactic template, supposed to cor- 
respond to a linear combination of synchrotron and free-free 
emission, is obtained from the difference between the K and 
Ka bands. This template is produced at one degree resolution. 
An additional free-free template is obtained from Ha emission 
(Finkbeiner 2003) corrected for dust extinction (Bennett et al. 
2003). A dust template is obtained from model 8 of Finkbeiner 
( 1999). "Clean" Q, V and W maps are obtained by decor- 



let al. 



relation of these templates from the original Q, V and W obser- 
vations. 

The main limitation of this approach is that the model used 
is insufficient to guarantee a good fit of the total foreground 
emission simultaneously inside and outside of the Kp2 mask. 
As a consequence, the maps produced are heavily contaminated 
by foregrounds in the galactic plane, the priority being given to 
higher galactic latitudes, with the objective of obtaining a part- 
sky high quality map on which high multipole CMB power spec- 
tra could be estimated reliably. 

In addition, the maps are likely to depend significantly on 
the prior model assumed. Here, the WMAP team chooses dust 
model number 8 of |Finkbeiner et al.|(|1999i, and also ignores 



the plausible existence of anomalous dust emission. The exact 
impact of these a priori decisions is difficult to evaluate. 

As an additional drawback, we note that the method gener- 
ates correlated noise in the foreground-reduced maps, originat- 
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Fig. 1. This figure shows a selection of presently available CMB maps from WMAP data. Contamination by galactic emission is 
visible in all of them at various levels, except for the 3-degree resolution EGS3 map (for which a galactic cut was applied and then 
filled by a plausible CMB extrapolated from higher galactic latitude data). 



ing either from K and Ka channel noise or from a background 
of radio sources. Finally, on supra degree scales, the K and Ka 
bands, which are the most sensitive ones, are used only to sub- 
tract foregrounds, whereas in the cleanest regions of the sky they 
would be more usefully used to estimate the CMB emission. 

For all these reasons, WMAP foreground-reduced maps are 
not good CMB maps according to the criteria listed in the intro- 
duction. 



An interesting remark from Eriksen et al. is that the amount 
of residual dust is high in the ILC maps - the method being able 
to subtract only about half of the dust present in the W band. 
At scales larger than 1 degree, this lack of performance is likely 
to be due to the part of dust emission uncorrelated to low fre- 
quency galactic foregrounds. On the smallest scales, the situa- 
tion is worse, as the low frequency WMAP channels do not have 
the resolution to help remove small scale dust emission from W 
band observations. 



3.1.3. The ILC by Eriksen et al. 



Eriksen et al. (2004) have obtained a CMB map at 1 degree res- 
olution with another implementation of the ILC. The map, de- 
noted here EILC1, uses only one year data. 



For this reason, in the present work, we improve on dust re- 
moval by using, as an additional measurement, the IRIS 100 mi- 
cron dust template obtained from a combination of DIRBE and 
IRAS maps ( Miville-Deschenes & Lagache 2005 i. As compared 
to the map of Eriksen et al., we also aim at better angular reso- 
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lution - and, obviously, better sensitivity, achieved by using five 
year data sets. 

3.1 .4. The Gibbs-sampling map by Eriksen et al. 



Recently, Eriksen et al. (2007) have produced a low resolution (3 
degree) CMB map from WMAP three year data, using a Gibbs 
sampling technique to explore the likelihood of a parametric 
model of CMB and foreground emissions. A CMB map is ob- 
tained as the average of 100 sample CMB maps drawn each at 
random according to the posterior distribution of the model pa- 
rameters given the observations. 

The free parameters in the model are spherical harmonic co- 
efficients acm of the CMB map, the CMB harmonic power spec- 
trum C[, monopole and dipole amplitudes in each WMAP band, 
the amplitude a(v) of a dust template in each band, and ampli- 
tudes f(p) and spectral indices /3(p) of a low-frequency fore- 
ground component, for each map pixel p. 

The model is constrained by fixing the dust template at 94 
GHz according to Fink beiner et al.| ([1999), by a prior on the 
low-frequency foreground spectral index, assumed to be close 
to that of synchrotron (-3 ± 0.3), and by the constraint that the 
monopole and dipole coefficients are orthogonal to the (noise- 
weighted) pixel-averaged foreground spectrum. 

In spite of a good fit of the assumed model to the data at 
high galactic latitudes, there are some strong limitations to the 
resulting CMB map, and hence to its usage: 

- The result of the sampling is obtained assuming a paramet- 
ric model of foreground emission. There is no possible way 
of validating the systematic errors due to mismodelling, ex- 
cept marginalising over all possible model skies. This would 
require a Monte-Carlo simulation which takes into account 
all uncertainties in the modelling, not only values of the pa- 
rameters for a given parametric model, but also the choice of 
the parameter set to be used to model the foregrounds (vary- 
ing the dust template according to uncertainties, assuming 
a different foreground model, etc.). This is not presently 
available; 

- The resulting map is at 3 degree resolution and HEALPix 
nside = 64, considerably worse than WMAP can do; 

- The data sets are cut with the Kp2 mask. Although a CMB 
is recovered in the mask by the average of the sample maps, 
the effective resolution inside the mask is lower than in the 
rest of the sky. In some sense, the Gibbs sampling technique 
(as implemented here) allows to recover in the mask what 
is predictable from the outside map (assuming stationarity 
of the CMB anisotropy field). It allows only for a limited 
prediction of the CMB signal in the masked zone. 

For all of these reasons, the Gibbs-sampling map of Eriksen 
et al. (hereafter EGS3) is not a good "best CMB" map according 
to our criteria. 



3.1.5. The ILC by Park etal. 

Park et al. ( 2007[ l provide their own version (hereafter the PILC3 
map) of a one degree resolution CMB map obtained by an ILC 
on WMAP three year data. The originality of their approach lies 
in the fact that they cut the sky into 400 pixel ensembles, selected 
from a prior on their spectral properties. The 400 ensembles are 
defined from 20 x 20 spectral index bins (20 for K-V spectral 
index, and 20 for V-W). This approach is motivated by the fact 
that ILC weights are expected to vary with varying foreground 
properties. 



There is a weak point to this approach. The authors use, 
to define their pixel 'bins', the MEM solution derived by the 
WMAP team. If the MEM solution is wrong for a given pixel, 
that pixel will automatically be classified in the wrong pixel en- 
semble, and be weighted using the weights of the wrong popula- 
tion of pixels. To some extent then, this binning forces the result 
of the ILC to match the prior assumptions given by the MEM re- 
sults. In turn, the MEM solution uses as a prior the result of the 
WMAP ILC, which is subtracted from the WMAP observations 
prior to using the MEM method to separate galactic foregrounds. 

As a consequence, the connection of the CMB map of Park 
et al. to the original WMAP data is far from direct. The map is 
bound to bear the signature of any arbitrary choice made before, 
in particular the choice of WMAP ILC regions, the 3-component 
model for galactic emission, and the MEM priors. For instance, 
discontinuities at the boundaries between the 12 regions of the 
WMAP ILC are clearly visible in the map of K-V spectral index 
used by the authors, as well as their group index (see Figures 3a 
and 4a of their paper). 

Park et al. then investigate the error in their reconstructed 
map by Monte-Carlo simulations. However, they use as an in- 
put galactic emission template the very model obtained by the 
MEM. This means that in the simulations, the spectral index 
maps are "exact". Therefore, the simulations investigate accu- 
rately the errors only if the MEM solution is correct, which is 
not likely to be the case -at least not to the level of precision 
required for producing a CMB map useful for precision cosmol- 
ogy. 

Our method described in Section |4] uses as little prior infor- 
mation as possible, and aims at better angular resolution than 1 
degree. 

3.1 .6. The "clean" map of Tegmark et al. 



The approach of Tegma rk et al.| ( |2003| l, on both one year and 
three year data, is the only work to date which aims at produc- 
ing a CMB map with both full sky coverage and best possible 
resolution. 

Teg mark et al.| ( |2003| l have performed a foreground analy- 
sis of the WMAP one year maps, producing two high resolu- 
tion maps of the CMB. The first one, the "clean" map (hereafter 
TILC1), is obtained by a variant of the ILC in which weights are 
allowed to vary as a function of {, and are computed indepen- 
dently in nine independent regions. The second, the "Wiener" 
map, is a Wiener-filtered version of the same map, in which the 
Wiener filter is applied in harmonic space, but independently in 
each zone. 

This work by Tegmark et al. is an early attempt at find- 
ing linear combinations of the WMAP data which vary both in 
harmonic space and in pixel space. The pixel variation of the 
weights is made using zones which are defined according to 
the level of contamination by foregrounds, as computed from 
WMAP map differences W-V, V-Q, Q-K and K-Ka. The authors 
do not specify exactly how the frequency range is divided into l- 
bands. Whereas the text seems to indicate that weights are com- 
puted independently for each E, figures hint that the weights are 
actually band-averaged, in 50 logarithmic bands subdividing the 
multipole range. In the end, the authors obtain a CMB "clean" 
map with a "beam corresponding to the highest-resolution map 
band", i.e. the beam of the W band. 

The original paper describes the work done on the one year 
WMAP data. However maps for three year data are available on 
Max Tegmark's web site (see Table [Til. We use the three year 
map (TILC3) for comparison with our own solution. 
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Although the approach of Tegmark and collaborators is quite 
good at high galactic latitudes, we can see on FigurefTlthat it per- 
forms poorly in the galactic plane. Also, the authors have not re- 
moved detected point sources from the WMAP data before mak- 
ing the ILC. As a result, their CMB map contains obvious point 
source residuals, for instance around galactic longitude 305° and 
latitude 57°, where a 5 mK peak can be seen. 

Our method, although bearing some similarity with that of 
Tegmark and collaborators, aims at improving significantly the 
error characterisation, as well as the quality of foreground clean- 
ing in the galactic plane. 



3.1 .7. The "Wiener" map of Tegmark et al. 



In addition to their TILC map, Tegmark et al. ( 2003 ) publish a 
Wiener map (hereafter TW map), obtained from the TILC map 
by independent Wiener filtering in the 9 regions. This results 
in reduced integrated error in all regions, at the price of pixel- 
dependent extra smoothing. The consequence of this filtering is 
an effective zone-dependent beam. 

Because of this extra smoothing, it is difficult to compare the 
TW map with other maps. The most meaningful figure of merit 
for the Wiener map, anyway, is the actual power of the error (out- 
put map minus true CMB). This is unavailable for any useful up- 
to-date real data set. Additional discussion about Wiener-filtered 
maps is deferred to Section [6] 

3.1 .8. The ILC map of Kim et al. 



More recently, |Kim et aL] ( |2"00 8) have made a CMB map from 
WMAP five year data, using an "harmonic" ILC method (KILC5 
hereafter). Their method perform an ILC in the pixel domain 
but with pixel-dependent weights. The ILC weights are not con- 
stant over predefined zones on the sky but are computed as 
smooth weight maps defined in terms of an harmonic decom- 
position (hence the qualification of the method). More specifi- 
cally, the weight maps are determined by minimising the total 
output CMB map variance with the constraint that these maps 
have no multipoles higher than Cutoff- For stability reasons, the 
KILC5 map is obtained with f cutoff = 7. Prior to computing ILC 
weights, all maps are deconvolved from their beam (effectively 
blowing up noise on small scales, in particular for the lowest fre- 
quency channels). Then, the channels are combined using map 
modes for I < 300. 

With the above choices, the reconstructed map can not be 
good on small scales. As the authors notice themselves, using 
small scale modes results in minimisation of noise rather than 
foregrounds (and, obviously, rejecting the low-frequency obser- 
vations, which are the noisiest on small scales after deconvolu- 
tion from the beam). Better results could probably be obtained 
by estimating weight maps for different bands of I. In essence, 
this is what our needlet ILC method permits to achieve. 

As a last comment, we note that limiting the number of 
modes of the weight maps to I < 7, reportedly for reasons of 
singularity of the system to be solved, results in spatial coher- 
ence of the weights on scales of about 35 degrees. The galactic 
ridge, however, is about 1 degree thick. Hence, the spatial vari- 
ability of the ILC weights achieved by |Kim et al.| ( [20 08) is not 
quite adapted to the actual scale of foreground variation. The 
needlet ILC method presented in our paper, as will be seen later 
on, solves this issue in a very natural way. 



3.1.9. Other maps 

Other authors have performed various foreground cleaning in 
the WMAP observations, producing a number of CMB maps 
for several different models of the foreground emission. Bonaldi 
et al. (2007) perform component separation on WMAP data us- 
ing the CCA method described in Bonaldi et al. (2006). Maino 



et al. (2007) obtain also several CMB maps, using the FastICA 



method (Main o et al.||2002[ l. None of the CMB maps obtained 
is full sky, nor publicly available yet. They are not considered 
further in this analysis. 

Finally, some foreground cleaning has also been performed 
by Saha et al. ( |2007| l. Their paper also includes an interesting 
analysis of the ILC bias. The primary goal of that work, however, 
is to compute the CMB power spectrum, rather than producing a 
CMB map. 

3.1.10. Existing map summary 

Figure [T] shows six available maps, all displayed in the same 
colour scale. It illustrates the resolution and foreground contam- 
ination of the various maps. Table [T] summarises the main prop- 
erties of the maps. Only the TILC1 and TILC3 maps are high 
resolution attempts at component separation everywhere, includ- 
ing the galactic plane, combining all WMAP observations. All 
other maps (WILC maps, EILC1, EGS3, PILC3, KILC5) are at 
reduced resolution. 



3.2. Map comparison 

CMB maps produced from WMAP one year data have been 
compared by Erik sen et al.| ( |2004| l, showing quite large differ- 
ences, ranging from -100 to 100/zK. Similarly, |Park et al. (2007) 
compare their PILC3 map with the WILC3 and TILC3 maps at 
1 .4 degree resolution, showing differences in excess of 40 /jK. 
In the following, discrepancies between these various solutions 
are further investigated. 

As a first step, we evaluate by how much the various maps 
at one degree resolution disagree on the full sky. The top panel 
of Figure[2]shows the pixel based standard deviation of five such 
maps. In this evaluation, we exclude the EGS3 map (which is at 
3 degree resolution, and is the source of additional variance in 
the Kp2 sky mask, where the CMB is estimated at even lower 
resolution). 

All maps are obtained from observation of the same sky 
(sometimes starting from the exact same data set) and are 
smoothed to the same resolution (one degree) where noise is 
small. Discrepancies originate essentially from systematic dif- 
ferences in the methods. In the central regions of the galactic 
plane, discrepancies are significantly higher than 50 fiK (with 
bright spots above 90 fiK). In other sky regions, they are typ- 
ically in the 10-20 fiK range, except on compact spots again, 
where they are above 50 fiK. The latter are probably due to resid- 
uals of the emission from strong compact sources. The discrep- 
ancy between the maps is larger in the ecliptic plane, which is 
a signature of the impact of the instrumental noise, uneven be- 
cause of the WMAP scanning strategy. This is due to the differ- 
ent attention given to minimising the contribution of instrumen- 
tal noise on small scales (rather than foregrounds) in the final 
map (and to a much lesser extent to the difference of noise level 
between the various WMAP releases). 

The bottom panel of Figure [2] focuses the comparison to 
larger scales (3 degree beam). At high galactic latitudes, the dis- 
crepancies are below 10 fjK except for a few localised regions 
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(LMC, Ophiuchus complex) where they reach about 30 fiK. 
Differences in the North Polar Spur, at the level of 10-15 fiK, 
are also clearly visible. Close to and inside the galactic plane, 
systematic discrepancies significantly exceed 30 fiK. 

Including the EGS3 map for comparisons at high galactic lat- 
itudes does not change these conclusions. Pairwise comparisons 
of the available maps typically show the same level of discrep- 
ancies, which indicates that the variance of the solutions is not 
due to one single map in strong disagreement with all the others. 







4. The ILC needlet method 

4.1. The choice of the ILC 

It is striking that all the presently available full sky CMB maps 
derived from an analysis of the WMAP data have been obtained 
by an implementation of the ILC method. 
The ILC, indeed, has many advantages: 

- The method relies only on two very safe assumptions: the 
CMB emission law, and the fact that the CMB template is 
not correlated to foreground emissionr] 

- Under these assumptions, the method minimises the empiri- 
cal variance of the reconstruction error; 

- The ILC is very simple in implementation. 

The ILC has also two major drawbacks. 

- As noted for instance by Hinshaw et al. ( 2007 ), Delabrouille 



& Cardoso (2007 1, Saha et al.| ( )2.007} >, empirical correlations 



between the CMB and the source of contamination results in 
a bias; this bias is discussed in more detail below; 
In absence of a model of the contaminants (foregrounds and 
noise), it is not possible to predict the reconstruction errors, 
which somewhat annihilates the benefit of making very safe 
assumptions about the properties of the data set. 



i m ■!■■ 




Fig. 2. Top: the standard deviation, per pixel for nside=512, of 
the EILC1, TILC3, PILC3, WILC5 and KILC5 maps, at 1° res- 
olution. Bottom: same at 3° resolution. Note the different color 
scales for the two panels. 



The conclusion of this comparison is that foreground resid- 
uals exist in the published CMB maps at the level of about 50 
to 100 fiK in the galactic ridge, 20 to 50 fiK at low galactic lat- 
itudes, and 10 to 20 fiK at higher galactic latitudes (above 30 
degrees). 

This observation calls for localised weightings, adapted to 
local properties of the foregrounds and the noise. This, however, 
is not easily compatible with the recovery of the largest modes 
of the CMB map, as pointed out before by |Eriksen et al. (2004), 
and as demonstrated by the discontinuities observed between the 
CMB solutions in the different regions when the sky is cut, as in 
the TILC, EILC1, and WILC maps. 

Our approach, then, will be to vary the relative weightings 
on small scale for small scale CMB reconstruction, and keep the 
weighting uniform over large regions of the sky for the recov- 
ery of the largest scales. This can be achieved quite straightfor- 
wardly by using the spherical needlets discussed in Section [4] 



4.2. The ILC bias 

The existence of a "bias" in maps obtained by an ILC method is 
a well established fact. The derivation of this bias (which is, in 
fact, the systematic cancelling of a fraction of the projection of 
the CMB map onto the vector space spanned by the noise reali- 
sations for all the considered input maps), is given in appendix. 

The order of magnitude to keep in mind is that about (m - 
1) "modes" of the original CMB, out of N p , are cancelled by 
the ILC, where m is the number of channels used, and N p the 
number of independent pixels or modes in the regions for which 
the ILC is implemented independently. Note however that when 
the signal is strongly correlated between pixels, the bias can be 
significantly larger -see the appendix for details. 

The practical consequences are: 

- a loss of CMB power, which has to be taken into account for 
power spectrum estimation; 

- an anti-correlation of the map reconstruction error with the 
real CMB sky. 

The level of the bias induced by our method is investigated 
both theoretically (Appendix[A]i, and through Monte-Carlo sim- 
ulations. 



4.3. Needlets 

A frame is a collection of functions with properties close to those 
of a basis. Tight frames share many properties with orthonormal 
bases, but are redundant (see Daubechies ( 1992) for details). 

Needlets were introduced by Narcowich et al. (2006) as a 
particular construction of a wavelet frame on the sphere. They 
have been studied in a statistical context (e.g. |Baldi et al.| ([2008 ); 
Baldi et al. (2007 )) and have also been used recently for cosmo- 
logical data analysis problems (e.g.|Pietrobon et al.|2006). The 



2 In reality, it is likely that the CMB map actually is somewhat cor- 
related to the foregrounds (extragalactic point sources and SZ effect), 
because of the ISW effect. The implication of this is not studied further 
in the present paper. 
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most distinctive property of the needlets is their simultaneous 
perfect localization in the spherical harmonic domain (actually 
they are spherical polynomials) and potentially excellent local- 
ization in the spatial domain. 

We recall here the definition and practical implementation of 
the needlet coe fficents, following the generalised formulation by 



Guilloux et al. 



%<.» 



(2008 ). Let h , j e J be a collection of window 



functions in trie multipole domain, indexed by j. Suppose that 
for each scale j, £r is a grid of points (indexed by k e K^) 

which satisfies an exacj 3 quadrature condition with weights X£\ 
The needlets are axisymmetric functions defined by 



^(0 



J>»I<<?- 



a. 



(6) 



f=o 



where L[ denote the Legendre polynomials. 

Any square integrable function / on the sphere can be anal- 
ysed by the scalar products pr := (/, ifrr) of the function / with 
analysis needlets. All the needlet coefficients of scale j are ad- 
vantageously computed in the spherical harmonic domain, as the 
evaluation at points ft 7 of a function whose multipole moment 

are simply h ( a( m . These needlet coefficients, denoted yP, are 
given by: 



7? 



W0? 



Each field of needlet coefficients can in turn be convolved 
with some synthesis needlets 



fiXd = 




^L e ^-^\ 



(7) 



using the exact same procedure and leading to the map X^ 
whose multipole moments are h ( J h ( J ae m . The analysis and 
synthesis operations are summed up as: 



Analysis: 



SHT 

A — > a e „ 



(;) SHI ( ■) 

hf a tm => fl' 



Synthesis: 



y? 



SHT" 



h) J 'a {m 



SHT 



X U) 



hfa tm -A hf^- 

Double arrows denote as many transforms as scales in J~. If X 
is band-limited to ( < f max and if the reconstruction condition 
Yjj nj hj - 1 holds for all £ < £ max , then the complete process 
yields a decomposition of X in smooth maps, namely 



V£ X©^X W © 



(8) 



Note that the existence of a fast inverse spherical harmonic 
transform using the quadrature points &' is required in prac- 
tice, and that HEALPix pixels and weights fulfil the quadra- 
ture condition only approximately. Further details can be found 
in |Guilloux et al.| {2008 ), with an extensive discussion on the 
choice of the spectral window functions. 

A key feature of the needlet decomposition follows from the 
localization of the analysis functions which allows for localised 
processing (such as denoising, signal enhancement, masking 
etc.) in the needlet coefficient domain, i.e. applying some non- 



uniform transforms to the coefficients 



y k 



.a) 



4.4. The method 

The method implemented in this work, and applied both to simu- 
lations and to the real WMAP data sets (for all releases), consists 
in the following steps: 

- We start with the data set consisting of band-averaged tem- 
perature maps from WMAP (simulated or real data), to 
which we add the IRIS 100 micron map; 

- WMAP-detected point sources are subtracted from the 
WMAP maps; 

- We apply a preprocessing mask, in which a very small 
number of very bright, compact regions, are blanked (see 
Table [2}; blanked regions are filled-in by interpolation; this 
is done only on the real WMAP data. 

- All maps are deconvolved to the same resolution (that of the 
W channeQ); this operation is performed in harmonic space; 

- Maps are analysed into a set of needlet coefficients yJ fol- 



lowing the method described in 4.3 



- For each scale, the covariance matrix R of the observations 
is computed locally (using an average of 32 x 32 needlet 
coefficients); 

- The ILC solution is implemented for each scale in local 
patches; 

- An output CMB map is reconstructed from the ILC filtered 
needlet coefficients; this map constitutes our main CMB 
product; 

- That map is Wiener-filtered in harmonic space, to make an 
alternate CMB map with lower integrated error (our best 
guess CMB map). 

- In parallel, the actual ILC filter used on the analysed data set 
is applied to 100 different simulations of the WMAP noise, 
to estimate the noise contribution to the final map. 

- The level of the biasing, which depends on the geometry and 
not much on the actual templates of CMB and foreground 
and noise, is estimated on a set of fully simulated data. 

Each of these steps is described in more detail below. 



4.5. Point source subtraction 

Strong point sources in the input data set typically leave de- 
tectable residuals in the output ILC map, and twist at the same 
time the estimation of the background, thus conducting to lower 
rejection of other contaminants. On the other hand, their specific 
shape usually allows effective estimation and removal by other 
methods. For the purpose of this study we used information from 
the WMAP source catalogue (Hinshaw et aL]( 2007[ l) which pro- 
vides characterisation for all point sources detected above a 5<x 
threshold away from the galactic plane. 

For all sources identified, we subtract from the input maps a 
Gaussian profile at the given position and with the given flux. 
Conversion factors between flux density and Gaussian ampli- 
tude, as well as the FWHM of the Gaussian profile are taken 
from Table 5 of Page et al. (2003) for one year data, and corre- 
sponding updates for the more recent releases. 

For the simulated data set, the subtraction of detected point 
sources is mimicked by removing from the simulation all sources 
above 1 Jy (independently in all channels). 



or almost exact, for all practical purposes 



4 A noise weighted average beam is obtained from the Wl W2 W3 
W4 beam coefficients provided by the WMAP team 
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Name 


Galactic coordinates 


Type 


Crab neb 


184.5575 -05.7843 


SNR 


sgr A 


000.064+00.147 


Radio-Source 


sgrB 


000.599 +00.002 


Radio-Source 


sgrC 


359.4288 -00.0898 


HII region 


sgr D 


001.131 -00.106 


Molecular cloud 


Orion A 


209.0137 -19.3816 


HII region 


Orion B 


206.5345 -16.3539 


Molecular cloud 


Omega neb 


015.051 -00.674 


HII region 


CenA 


309.5159+19.4173 


QSO 


CasA 


111.735-02.130 


SNR 


Carina neb 


287.6099 -00.8542 


HII region 



Table 2. List of compact regions blanked in the pre-processing 
step. A circular patch centered on the source, of radius 75, 55, 
45, 45, 34 and 35 arcminute for the K, Ka, Q, V, W and IRIS 100 
ymi bands respectively, is masked. The masked regions are then 
filled with an extrapolation of edge values. 



Band (J) 


tmin 


Mnax 


nside(j') 


1 





15 


8 


2 


9 


31 


16 


3 


17 


63 


32 


4 


33 


127 


64 


5 


65 


255 


128 


6 


129 


511 


256 


7 


257 


767 


512 


8 


513 


1023 


512 


9 


769 


1199 


512 



Table 3. Spectral bands used for the needlet decomposition in 
this analysis. Needlet coefficient maps are made at different val- 
ues of nside, given in the last column. 



power of 2 larger than / max /2, with a maximum of 512. Details 
about the bands used are given in table [3] and figure [3] 



4.6. Blanking of compact regions 

In addition to the point sources subtracted above, some compact 
regions of strong emission (mostly in the galactic plane) exceed 
the rejection capabilities of the method used in this analysis, be- 
cause they are too local and/or too specific. Their contribution in 
the wings of the needlets also contaminate the solution far from 
the centre of the sources. Those sources cannot be satisfactorily 
subtracted in the same way than the previous ones, either be- 
cause they are not strictly speaking point-like, or because they 
are bright enough that small departures of actual beam shapes 
from the Gaussian model used in the subtraction step leave sig- 
nificant residuals. As they represent only a very tiny fraction 
of the sky (we single out eleven such sources), we blank out 
these regions in all WMAP channels, cutting out circular patches 
adapted to the size of the beam and of the source. Table P2] gives 
the list of those regions with their main characteristics. 

To reduce local pollution of the needlet coefficients by the 
sharp cut, the small blanked regions are filled in by a smooth in- 
terpolation, so that fluctuations at a larger scale than the hole size 
are at least coarsely reconstructed. More precisely, interpolation 
is made by diffusion of the boundary values inside the hole. 

Although this masking and interpolation has no reason to 
be optimal, it is an efficient way of reducing the impact of 
very strong sources on their environment. The CMB inside the 
masked patches is recovered (to some extent) both by the inter- 
polation of original maps (which avoids sharp discontinuities) 
and by the needlet decomposition and ILC reconstruction. The 
masked region is tiny: 0.058% of the sky in the K channel (the 
most affected). 



4. 7. Needlet decomposition 

The original observations (WMAP and IRIS) are decomposed 
into a set of filtered maps represented by their spherical har- 
monic coefficients: 
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where at, 



h { ti{ m 



(9) 



are the spherical harmonics coefficients of the origi- 
nal map, and aP those of the same map filtered by the window 

function j. Needlet coefficients y k J are obtained as the value of 
the filtered map at points &. 

For each scale /, the coefficients yr are computed on a 
HEALPix grid at some value of nside, compatible with the 
maximum value of 6 of band j. We use for nside(j) the smallest 
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Fig. 3. The spectral bands used in this work for the definition of 
the needlets. 



4.8. ILC implementation on needlet coefficients 

The general idea is to implement independently the ILC on sub- 
sets of the needlet coefficients yr . For a given scale, these coef- 
ficients come in the format of a set of HEALPix maps (one per 
frequency channel). The ILC is implemented locally in space 
and locally in I as follows. 



j(/3 



U)AJV 



Covariance matrices Re; = (yr yr ) for scale j at pixel 
k are estimated as the average of the product of the computed 
needlet coefficients over some space domain D^. In practice, for 
the present analysis, we make use of the hierarchical property of 
the HEALPix pixellisation, and compute covariance matrices as 
the average in larger pixels, corresponding to a HEALPix pixelli- 
sation with nside = nside(/)/32. This provides a computation 
of the statistics by averaging 32 2 = 1024 samples, which results 
in a precision of order 3 per cent for all entries of the R^ matrix. 
It implies an ILC bias of order 5/1024 (for m = 6 channels and 
N p = 1024 coefficients per domain on which the ILC filter is 
estimated independently (see appendix for details). 

We denote as R^, the estimate of RJr obtained by averag- 
ing the value of yr yr in domain D^. 

On the largest scales (I < 50), the typical angular extent of a 
needlet is larger than 5 degrees, and the value of nside for the 
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map of needlet coefficients is less than 32. Covariance matrices 
are then computed on the full sky rather than on the largest pos- 
sible HEALPix grid, i.e. Dt is the complete sky, rather than one 
of the 12 basis Healpix pixels. 

Using these covariance matrices, the ILC is implemented us- 
ing equationfflfor each domain. The estimated CMB needlet co- 
efficients are: 



R1 



CMB 



* T HU 



■7? 



(10) 



4.9. Full map reconstruction 

The full CMB map reconstructed from this set of needlet coeffi- 
cients is our basic needlet ILC (NILC) CMB map. 



4.10. Final Wiener filtering 

For a number of purposes, in particular subtraction of an esti- 
mate of the CMB to study other emissions, it is interesting to 
use, instead of our ILC map at the resolution of the WMAP W 
channel, a map with minimal error. Such a map is obtained from 
the ILC map by one-dimensional Wiener filtering. 

As a last processing step towards a minimum variance CMB 
map, we thus Wiener-filter our CMB map, to get rid of the large 
noise contamination at high I. The Wiener filter is performed in 



i.e. W{ — 



harmonic space as described in 2.3.2 

The harmonic Wiener filter is given by formula [5] 
b 2 ( C{l{b 2 ( C{ + N(). For its implementation, we need to know the 
relative power of CMB and noise. We assume that the best fit 
CMB power spectrum of the WMAP team is correct, hence C( 
is known. The beam factor b c is assumed perfectly known as 
well. The denominator b 2 Ci + N( can be estimated directly as 
the power spectrum of our output needlet ILC map. 

In practice, we smooth the power spectra with 6£/£ — 0. 1 to 
lower the variance of the power spectrum estimator on the output 
needlet ILC map. Even with this, the filter is poorly estimated for 
low modes, because of the large cosmic variance. As can be seen 
on Figure|5] the signal to noise ratio of our reconstructed map is 
expected to be quite high at low I. Therefore, the Wiener filter 
for low modes is expected to be very close to 1 . For this reason, 
we set W( - \ for t, < 200, and use a linear interpolation between 
I = 200 and £ = 250. 



4.11. Noise level estimate 

The level of noise contamination (variance per pixel, and average 
power spectrum) in the output map is estimated by Monte-Carlo, 
using the average of 100 realisations of the WMAP noise maps. 
For each initial set (/) of five noise maps (one noise map per 
WMAP channel), a single output noise map rip is obtained by 
performing the needlet decomposition of the initial noise maps, 
and filtering needlet coefficient maps with the same filter as that 
used on the single full simulated data set. 

Denoting as n p and ny respectively the pixel value and the 
harmonic space value of the noise map number i, we compute: 
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These are respectively estimates of the noise pixel variance 
and of the noise power spectrum of our final map. 

4.12. Bias estimates 

The impact of the ILC bias is estimated by Monte-Carlo simula- 
tions on simulated data sets. The corresponding error is of order 
2% of the CMB. 



4.13. Practical implementation 

The practical implementation of this processing pipeline is made 
essentially using the octave language (the free software version 
of Matlab). The analysis is done in the framework of the pipeline 
tool developed in the ADAMIS team at the APC laboratory. This 
tool provides a flexible and convenient web interface for running 
our data analysis on simulations or real data with easy handling 
and tracing of the various pipeline optionslj Single runs of the 
full pipeline require less than half an hour on a single processor 
of a standard desktop computer (dominated by harmonic trans- 
forms), whereas numerous pipelines on simulated data sets for 
Monte-Carlo are run on the ADAMIS 88-processor cluster, op- 
timised for efficient I/O. 



5. Simulations 

5. 1 . The simulated data 

We start with a validation of our method on simulated data sets. 
For this experiment, synthetic observations of the sky emission 
are generated using the Planck Sky Model (PSM). The PSM is 
a flexible software library, designed for simulating the total sky 
emission in the 10-1000 GHz frequency range, and developed as 
part of the foreground modelling activities of the Planck work- 
ing group on component separation (Planck WG2). Sky emis- 
sion comprises galactic components of four origin (free-free, 
synchrotron, thermal dust, and spinning dust, with spectral emis- 
sion laws for dust and synchrotron varying from pixel to pixel), 
CMB, kinetic and thermal SZ effect, and the emission from a 
population of galactic and extragalactic point sources which in- 
cludes radio sources, infrared sources, and an infrared back- 
ground. Although not perfect, the model sky is thought to be 
sufficiently representative of the complexity of the real sky emis- 
sion for our simulations to be meaningful. 

Sky maps are produced at WMAP central frequencies for 
the K, Ka, Q, V and W band, and convolved in harmonic space 
with approximate WMAP instrumental beams (Gaussian sym- 
metric beams are used for these simulations). Uncorrelated, non- 
stationary Gaussian noise is added, with a pixel variance de- 
duced from the WMAP sensitivity per channel and effective hit 
count. To mimic the subtraction of the brightest point sources 
detected by WMAP, we remove from the model sky, at each fre- 
quency, all sources with flux above 1 Jy (assuming they would 
have been detected, and can be subtracted from the data set). The 
1 1 compact regions listed in table [2] however, being specific to 
the real sky, are not blanked for the simulations. 

Although these simulations provide only an approximation 
of the real WMAP data sets, they are representative enough that 
the simulated data offer a component separation challenge close 
to that of the real data set. The IRIS map is used as part of the 
full set of data for the ILC implementation on simulations. 



see http://www.apc.univ-paris7.fr/APC_CS/Recherche/Adamis/ 



the 'outreach' section 
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power spectrum (dotted line), the spectrum of the output CMB 
(solid black line), and the spectrum of the map of residuals (dif- 
ference between output and input, dashed line), both full sky (top 
panel) and in the HGL region (bottom panel). The angular power 
spectrum of the residual map is seen to be small as compared to 
the CMB power on large scales, the two being comparable at 
I - 500. Noise dominates on smaller scales. The residuals due 
to the presence of galactic emission are seen to contribute power 
essentially below t - 400, where the power of the difference 
map is seen to be slightly higher in the full sky power spectrum 
than in the HGL power spectrum (this is visible, in particular, at 
the top of the first acoustic peak). 



3.02b 




\ rii 




Fig. 4. Top: the simulated input CMB map. Middle: the recon- 
structed CMB. Bottom: the difference (output-input), displaying 
the residuals left by the method. All three maps share the same 
colour scale, and are at the resolution of the WMAP W channel. 



5.2. Results 

Figure [4] shows the input simulated CMB, the output CMB, 
and the difference of the two for one particular simulation. 
The reconstruction is visually good except in regions of local 
strong galactic emission (in the galactic ridge, for example). This 
is to be expected: not only the method can not remove fore- 
grounds perfectly, but in addition the price to pay to remove 
foregrounds (even imperfectly) is more noise (because of sub- 
optimal weighting of the observations as far as noise contamina- 
tion is concerned). 

A more quantitative estimate of the level of contamination 
of the CMB map by foregrounds and noise is obtained by look- 
ing at power spectra. Figure [5] shows the input simulated CMB 
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Fig. 5. Top: For simulated data sets, full sky power spectra of 
output CMB map (plain line), and of difference map (dashed 
line). The CMB model used in the simulation is over-plotted as 
a dotted line. Bottom: Same at high galactic latitude only (HGL 
region). 



5.3. Bias 

As discussed in |4.2| we expect a (small) bias in the ILC map, 
due to empirical correlations between the CMB emission and 
contaminants (including noise and foregrounds). This is not par- 
ticular to our approach, and is expected for any ILC method. For 
better characterisation of our output map, we evaluate the effect 
both theoretically (in appendix) and numerically. 
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Although a general analytic estimate of the bias is compli- 
cated, appendix |A| shows that (to first order at least and for rea- 
sonable assumptions about the CMB, foregrounds, and noise) 
the amplitude of the effect does not depend much on what the 
actual foregrounds are in detail, but is set essentially by the ge- 
ometry of the domains considered (through a number of effective 
modes). It is then possible to estimate the amplitude of the effect 
by Monte-Carlo simulations on synthetic data sets resembling 
the actual WMAP observations. 

Figure [6]illustrates an estimate of the bias b(() as a function 
of the harmonic mode, computed as a fractional error: 
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where ai m are the harmonic modes of the input CMB map, and 
~5{ m the harmonic modes of the output CMB map. The numera- 
tor in this equation computes the covariance of the residual map 
and the input map as a function of (, and the denominator is a 
normalisation factor. For an error uncorrelated with the input, 
b{€) should be close to on average. Analytical estimates of the 
effect (see appendix) suggest a bias of order 2% for our imple- 
mentation. The numerical estimate of Figure [6] obtained as the 
average bias for 500 simulated data sets, is in good agreement 
with this prediction, with slight variations due to varying num- 
ber of effective modes for different needlet scales. The average 
value of b(€) between 1-2 and 1000 on that simulation is about 
-2.2 %. 




Fig. 6. The fractional bias in the ILC map as a function of (, 
for one single simulation. This figure is illustrative both of the 
amplitude of the effect, and on its variance for one single reali- 
sation. Bias, and standard deviation of the bias, are of the same 
order of magnitude for most of the I range. 



6. Application to WMAP data 

6.1. ILC Result 

We now turn to the description of the results obtained on the real 
WMAP data sets. In order to facilitate the comparison with exist- 
ing maps, we process independently one year, three year and five 
year data, to obtain three CMB maps (hereafter NILC1, NILC3 
and NILC5). For each year, we use the beam estimates, noise 
level, and point source catalogue provided with the correspond- 
ing release. 



The improvement of CMB reconstruction with consecutive 
data releases is illustrated on Figure IT] which shows the full 
sky power spectra of the NILC CMB maps obtained with one 
year, three year, and five year WMAP data. The power spectra 
displayed are the raw power spectra of the output map, com- 
puted full sky, and smoothed with a variable window in £ of 10% 
width. Whereas the lower part of the spectrum, cosmic variance 
limited and CMB dominated, does not change much, the high 
I spectrum of the map, dominated by noise, decreases substan- 
tially with increasing observation time - as expected. The ex- 
cellent agreement at low I (up to I =* 300) between the power 
spectra and the model is striking. The bumps in the spectrum due 
to the first and second acoustic peaks are cloearly visible on the 
five year map spectrum. 
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Fig. 7. Power spectra of the NILC map for one year, three year, 
and five year WMAP data 



Our full-sky reconstructed CMB map for the five year obser- 
vations, at the resolution of the W-channel, is displayed in the 
top panel of Figure [8] 



7. Discussion 

7.1. Comparison with other maps 

A full comparison of our needlet ILC maps (for all data re- 
leases) with all the other available maps would be too long for 
the present discussion. Rather, we decide to compare our five 
year result only with the TILC3 map on small scales (choice is 
motivated by the fact that the TILC is the only other full sky 
high resolution map available), and with the EGS3 map on large 
scales. This is of particularly interest, as the EGS3 is the only 
map obtained with a method not based on the ILC, and also is a 
method specifically implemented for the recovery of the largest 
scales. 



7.1 .1 . Comparison at the pixel level - small scales 

On the smallest scales, we compare our needlet ILC map with 
the TILC and with the WMAP foreground-reduced W band. 
Figure [10] shows local regions of the foreground-reduced map, 
the NILC5 map, and the TILC3 map, in the galactic plane and at 
the galactic pole. Our needlet map is clearly less contaminated 
by galactic emission than the other two. At high galactic lati- 
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Fig. 8. Top: the NILC5 reconstructed WMAP CMB at the reso- 
lution of the W channel. Middle: the harmonic Wiener NILC5 
CMB map. Bottom: The NILC5 CMB map at 3 degree resolu- 
tion. 



tude, the NILC5 and TILC3 are visually comparable, while the 
foreground-reduced map appears to be more noisy, as expected. 
The power spectrum of the output map for the five year data 
(NILC5 map), for three different sky coverages, is shown on the 
bottom panel of Figure [9] On the same panel, we plot the an- 
gular power spectrum C( corresponding to the WMAP best fit 
model, corrected for the W-channel beam. On the top panel of 
the same figure, we show the same power spectrum estimates 
for the TILC3 map. This shows the improvement achieved by 
our method close to the galactic plane. This improvement is due 
both to the needlet approach and to the use of the IRIS map to 
help with dust subtraction. As seen in Figure 17] the difference in 
quality between NILC5 and TILC3 cannot be explained solely 
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Fig. 9. Power spectrum of the reconstructed WMAP CMB map. 
For each of the two panels, the CMB best fit model is shown as 
a solid black line, and power spectra computed at low galactic 
latitudes (using the LGL mask), on the full sky (no mask), and at 
high galactic latitudes only (HGL mask) are displayed as dashed 
lines. Note that the map spectra plotted here are directly those 
of the maps, without any correction for the beam. Top panel: 
ILC map of |Tegmark eTal] ( [2003] ) (TILC3). Bottom panel: this 
work, with five year WMAP data. The power spectrum of the 
needlet ILC CMB is significantly more homogeneous than the 
power spectrum of the TILC3 map. We interpret this difference 
as an indication that the TILC3 map is significantly more con- 
taminated by residuals of galactic emission. Note the different 
scales of the y-axis, and the improvement on small scales, with 
a noise power of about 0.024 mK 2 at t = 1000 for the NILC5 
map at high galactic latitude, instead of about 0.040 mK 2 for 
TILC3. As indicated by Figure [JJ this is due essentially to the 
better quality of the five year release, since the NILC3 map also 
has a noise power spectrum of about 0.040 mK 2 at £ - 1000. 



by reduced noise (NILC5 and NILC3 being very close in quality 
for all scales except the smallest ones). 

7.1 .2. Comparison at the pixel level - large scales 

Figure [TT| gives a visual comparison of NILC5 (this work) and 
the EGS3 priksen et al.|2007] l, as well as of NILC5 and KILC5 
(Kim et al. 2008). On the top row, we display the EGS3 and 



KILC5 at a resolution of 3 degrees, and degraded to nside= 64. 
The bottom row shows the difference between our needlet ILC 
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Fig. 10. CMB maps from WMAP three year data obtained with three techniques: left column: WMAP foreground reduced, W 
channel; middle column: our needlet ILC map; right column: the TILC map. The top line corresponds to a patch located in the 
galactic plane, centred around coordinates (/, b) = (45°, 0°). The bottom line shows the recovered CMB around the North Galactic 
Pole. 



solution (displayed on the bottom panel of Figure [Sj and these 
two maps. 

The most striking difference between the five year needlet 
ILC map and the EGS3 is in the galactic plane, where the EGS3 
does not recover the intermediate angular scales. The difference, 
however, shows no clear particular structure, as expected if it 
is the random realisation of a Gaussian random field. It is thus 
probably essentially due to the difference between our CMB re- 
construction on scales larger than 3 degrees, and the CMB on 
larger scales that can be inferred from the CMB reconstructed 
by Eriksen et al. outside of their galactic mask. At higher galac- 
tic latitude, the two maps are in good agreement, with no obvious 
feature which could be correlated to foregrounds or to the CMB 
itself, with the exception of a hot spot in the large Magellanic 
cloud (which might be residual emission of the LMC in our map, 
as Erikse n~et al.| (2007) actually mask the centre of the LMC 
and obtain a solution in the direction of the LMC by extrapo- 
lation from nearby pixels). Above 30 degree absolute galactic 



latitude, the RMS of the difference between our 3 degree map 
and the EGS3 map is 5.7 //K. The two maps are in much better 
agreement than the EGS3 and the WMAP MEM model maps 
(see figure 3 in Eriksen et al.). Note however that theoretically, 
if there were no foregrounds in the WMAP data, the noise stan- 
dard deviation cr„ on a 3 degree map obtained by noise-weighted 
averaging using all WMAP channels would be about 3.2 jjK for 
three year data. If instead we assumed that only the three high- 
est frequency channels are free of foreground contamination, <r n 
would be 4.4 /J.K. 

The difference between our map and the KILC5 map is more 
systematic with, in particular, stronger differences in the galactic 
plane, in spite of the fact that the two methods work on the ex- 
act same input data set. A careful visual inspection of the CMB 
maps themselves gives the impression that the KILC5 map is 
probably systematically negative towards the galactic central re- 
gions. There is, however, no secure way to be certain which map 
is best. 
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Fig. 11. Top left: The EGS3 map. Top right: The KILC5 map, smoothed at 3 degree resolution. Bottom left: The difference NILC5- 
EGS3; we see is a clear difference in the galactic plane with no particular structure, compatible with a smooth Gaussian field, where 
the EGS3 solution poorly estimates intermediate scales; a patch in the difference map at the location of the Large Magellanic Cloud 
is clearly visible also. Bottom left: The difference NILC5-KILC5; a clear structure aligned with the galactic plane is clearly visible, 
with a big difference towards the galactic center. 



7.2. Map characterisation 

7.2.1. Beam 

The effective beam of the reconstructed maps are plotted in 



Figure 12 for both the full resolution five year needlet ILC map, 
and for the Wiener-filtered version. The map has been recon- 
structed for the range < I < 1200, with a smooth transition of 
the response, between I of 1024 and 1200, from the nominal W 
band beam value to 0. This smooth transition permits to avoid 
ringing effects which happen in case of a sharp cut-off in 6. The 
ratio of the Wiener beam (dashed line) and ILC beam (solid line) 
gives a measure of the signal to noise ratio in each mode. 

7.2.2. Instrumental noise 

Given the ILC filter computed on the real data set, the level 
and properties of the instrumental noise can be straightforwardly 
computed by applying the same filter to simulated WMAP noise 
maps. From 100 noise realisations, we compute the average full- 
sky noise power spectrum (Fig ure [L3) l, as well as the noise stan- 
dard deviation per pixel (figure |14[). Noise properties are not as 
simple as one may wish: the noise is non stationary, because of 
both the uneven sky coverage and of the localised processing. 
It is also somewhat correlated pixel-to-pixel, in particular close 
to the galactic plane. This is unavoidable, but fortunately our 
pipeline permits to make as many Monte-Carlo realisations of 
the instrumental noise as needed for any scientific study made 
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Fig. 12. Harmonic response of the beam of the NILC5 CMB map 
and of the Wiener-filtered version. 



using our needlet ILC map. Hundred such simulations are made 
available as part of our main data products. 



7.2.3. Foreground residuals 

More problematic is the evaluation of the contamination of the 
CMB map by foreground residuals. It requires prior informa- 
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Fig. 13. Plot of the power spectrum of the noise (solid line). The 
spectrum for the WMAP best fit model is shown as a sold line 
also, for comparison. The dashed line is 2% of the WMAP best 
fit C(, indicative of the level of the expected ILC bias. The bias is 
seen to dominate on large scales. There is, however, little margin 
for improvement, as few independent modes (or needlet coeffi- 
cients) are available on the largest scales. 




Fig. 14. Map of the standard deviation of the noise, per pixel at 
nside=512, for the five year needlet ILC map at the resolution 
of the W channel. The adaptation of the filter to local contamina- 
tion is obvious from the uniform noise level in large healpix pix- 
els, in particular in the galactic plane. These large healpix pixels 
correspond to nside=16, which is the size used for computing 
the average filter for the smallest scales. This illustrates the im- 
pact of the compromise between subtracting foregrounds and in- 
creasing the noise with the use of the lower frequency channels, 
which have poor sensitivity on small scales. 



tion about the foregrounds, which the ILC method avoids using. 
An indication of the level of systematic uncertainty is obtained 
from the comparison of the various solutions (figure[2]). An other 
option consists in checking the contamination on simulations. 
Figures |4]and[5]give an idea of the expected contamination from 
such simulations. This, however, is good only as long as the sim- 
ulations are representative enough. 

The comparison of the power spectrum of the output CMB 
map computed at high and low galactic latitude (figure |5J, and 
a visualisation of the output CMP map at high and low galactic 
latitude (figure 10 1 also give an idea of the amount of galac- 
tic residuals, but none of these estimates is fully satisfactory for 



careful CMB science. This is, however, not particular to our map. 
No published CMB map is provided yet with a good estimate 
of foreground contamination. Although this, clearly, is not fully 
satisfactory, we leave further investigations on this question for 
future work. 



7.3. Final comments 

7.3.1 . Is our map optimal? 

In the present work, we have obtained a CMB map which has 
been shown to be significantly less contaminated by foregrounds 
and noise than the other existing maps obtained from WMAP 
data. A natural question is whether we can do even better. 

In the following, we outline where there is margin for im- 
provement, and explain why we have stopped the analysis at its 
present state. 

First of all, the present analysis uses only limited external 
information and data sets: WMAP point source detections, and 
the IRIS 100 micron map. It is likely that something could be 
gained by using additional observations to help constraining the 
galactic emission. 

Second, there is a trade-off between localisation of the es- 
timation of covariance matrices, and bias in the ILC. The esti- 
mation of covariance matrices R v over sets of N p - 32 x 32 
needlet coefficents results from a compromise which has been 
made based on some attempts with varying N p on simulations, 
but has not been optimised in any way. In addition, the optimal 
solution is probably different at high galactic latitudes, where 
weights given to different channels probably do not vary much 
and thus require less localisation, and at low galactic latitudes, 
where the complexity of galactic emission calls for more locali- 
sation. We have tried to use N p - 16 x 16 (more localised ILC 
filters, but more bias) and N p = 64x64 (less localised ILC filters, 
but less bias). Our choice of N p = 32x32 seems, on simulations, 
not worse than anything else (nor much better either). The bias 
error has been verified to remain below the reconstruction error 
due to the contribution of the noise for small scales, and below 



cosmic variance uncertainty for large scales (see figure 13 I, and 
remains below an acceptable level of few percent. 

Similarly, the choice of the spectral window functions used 
on this data set has not been the object of specific optimisation. 
At low (, we follow a "dyadic'"scheme, where each window 
reaches an t m . dx about twice the previous one. Wide spectral win- 
dow functions allow for more localisation in pixel space, but nar- 
row window functions allow for more accuracy in the harmonic 
domain. At high I, because of the variation of the beams with 
t, the relative noise levels of the different channels change quite 
fast with I, which calls for more localisation in I space. Here 
again, the optimal window functions are probably not the same 
at high and at low galactic latitudes. In practice, we chose a small 
number of bands to limit the number of harmonic transforms in 
the pipeline and allow reasonable localisation of the analysis. 

The choice we have made results from the principle of sim- 
plicity. We have tried to devise a pipeline which depends as little 
as possible on external data, on priors, or on fine tuning. A very 
simple scheme has permitted to obtain a CMB map convincingly 
better than other maps available. This does not preclude any at- 
tempt at more optimisation for future work if needed. 

7.3.2. Why ILC and not ICA? 

It is certainly possible to tune our pipeline, changing some of its 
parameters. It would be possible also to use other methods than 
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an ILC, for instance Independent Component Analysis (ICA) 
metho ds such as SMICA ( [Delabrouille et al.|2003)|Cardoso et al.| 
2008), or more generally maximum likelihood methods fitting 
parametric models of the foregrounds on the data sets. Such 
methods extract information about the foregrounds from the data 
directly, possibly with the help of ancillary data sets, and use this 
information to clean the observations in some way. 

In the present case however, it is not very likely that such 
attempts would give much better results than what is obtained 
here, unless one uses a very significant number of additional 
data sets and safe prior information. Indeed, the WMAP data 
consists in five channels only, from which a component separa- 
tion method based on a meaningful model of foreground emis- 
sion needs to extract CMB, synchrotron, free-free, spinning dust, 
thermal dust (i.e. five templates), and possibly also point sources, 
and variations of the spectral indices of some of the components 
- not to mention special regions of galactic emission as cold 
cores and H-II regions, nor particular objects as nearby resolved 
galaxies or galaxy clusters. Any component separation method 
based on the estimation of parameters for such a rich model 
would almost certainly be confronted to indeterminacy issues. 
Methods based on a precise model are expected to be effective 
when the data are very redundant as compared with the num- 
ber of "parameters" of the emission model, which would not be 
the case here. Hence, the ILC is probably one of the best ap- 
proaches for doing component separation on WMAP data. It is 
not surprising, then, that all methods producing full sky CMB 
maps from WMAP, or nearly so, are implementing some vari- 
ant of the ILC. Incidentally, ICA methods could benefit from the 
needlet framework. 



8. Conclusions 

In this paper, we have described a new approach for implement- 
ing CMB extraction in WMAP data, using the ILC method on a 
needlet frame. Tests on simulations show excellent performance 
of the method, thanks to localisation both in pixel space and 
in harmonic space. Localisation in pixel space allows the ILC 
weights to adapt themselves to local conditions of foreground 
contamination and noise, whereas localisation in harmonic space 
allows to favour foreground rejection on large scales (where 
foregrounds dominate the total error) and noise rejection on 
small scales (where foregrounds are negligible but where, after 
beam deconvolution, the relative noise level between the various 
WMAP channels varies a lot as a function of scale). Needlets 
permit to vary the weights smoothly on large scales, and rapidly 
on small scales, which is not possible by cutting the sky in zones 
prior to any processing. 

As a further improvement to previous work on WMAP data, 
we include a dust template in the set of analysed observations. 
This is motivated by the fact that on the smallest scales, observed 
with reasonable signal to noise ratio by the W channel only, dust 
emission contributes a significant fraction of the total reconstruc- 
tion of the map. Using the IRIS 100 micron map as an additional 
observation enables the ILC to reduce the final contamination 
by dust -thanks to correlations of dust emission between the W 
channel and the 100 micron map. Special care was also taken 
to subtract from the data, prior to making the ILC, a number of 
strong point sources. 

As discussed at length in the main text and in the appendix, 
the implementation of a filter (the ILC) based on empirical esti- 
mates of covariance matrices leads to a bias. This is not partic- 
ular to our map, but is the case for any ILC map. We have es- 
timated the level of this systematic effect, both analytically and 



numerically, to be at the level of a few per cent on all scales. Our 
simulation tool permits to make accurate estimates of the ampli- 
tude of the effect, if needed for any scientific exploitation of the 
NILC5 map. 

The application of the method to WMAP one year and three 
year data (in addition to five year data) permits to compare the 
needlet ILC solution to previous work. Our map is seen to be at 
least as good as others on large scales, while being significantly 
less contaminated by residual foregrounds and noise than others 
on small scales, in particular in the vicinity of the galactic plane. 
The application of the method to WMAP five year data yields a 
CMB map which we believe to be the cleanest full sky map of 
the CMB to date. Contamination by noise, and the power loss 
due to the use of the ILC method, are characterised by means of 
Monte Carlo simulations. The map is available for download on 
the AD AMIS web site jj and can be used for a variety of science 
projects relying on accurate maps of the CMB. 
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Appendix A: Derivation of the ILC bias 

In this appendix, we compute the error made after CMB recon- 
struction with the ILC and some of its statistical properties. In 
particular, we derive the correlation of the error with the true 
CMB signal, which yields a non unit effective "response" of the 
ILC filter - and hence a bias in the reconstructed map and in the 
CMB power spectrum computed from it. This bias has to be ac- 
counted for properly for further uses of the reconstructed CMB 
map. 

We model the data as: 

x p - as p + n p (A.l) 

where s p are the coefficients of the map of interest over some 
domain T> (e.g. needlet coefficients of the CMB map for a given 
scale and a given patch of the sky, or pixel values in a certain 
region of the sky, or values of the harmonic coefficients in some 
band of I), p indexes the coefficient (i.e. pixel, or harmonic 
mode, or needlet coefficient). x p are the observations of coef- 
ficient p for the set of available observed maps, and n p the cor- 
responding "noise" terms (including foreground contaminants). 

The ILC is best applied over domains of p where all coeffi- 
cient have (near) uniform expected signal and noise properties, 
so that the ILC weights are optimal simultaneously for all p. In 
particular, the RMS values of all maps do not depend (much) on 
p in a given domain. Hence, we concentrate on one given domain 
of p for which we assume that the sequences s p and n p are inde- 
pendent, Gaussian random variables with distribution 7V(0, cr 2 ) 
and W(0, R„), with cr 2 the variance of s p (the CMB) and R„ the 
co variance matrix of the noise (including foregrounds). 

The ILC estimate of s p in domain T> is given, for all p, by 



„ a 1 R- 1 

s p ~ ,-=■ , x p 
a'R~ l a 



(A.2) 



where R x is an estimate of the covariance matrix of the observa- 
tions, obtained as: 



K v — 7 XnX„ 

N Lu p p 

' p 

= jT- ^(as p + iip)(as p + n p y 



(A3) 



with N p the number of coefficients in domain D. In the limit 
of large N p , R x approaches its expectation (ensemble average) 
value E(R t ) = R x , For finite N p , we have instead 



R. v = R x + A A . 



(A.4) 



where A x is a correction term corresponding to the departure of 
the empirical correlation from its ensemble average due to the 
finite sample size N p . From now on, we assume that N p is large 
enough that this correction is small: we investigate effects at first 
order in 1/N P . 

A.1. First order expansion of the reconstruction error 
We are interested in the reconstruction error: 



dp — Sp Sp 



(A.5) 



The first and second moments of d p , i.e. the mean value E(d p ) 
of the reconstruction error, as well as its variance E(d p 2 ), are of 



particular interest for the interpretation of the reconstructed map. 
In particular, we have 



E(^) = E(s p 2 ) + E(d p 2 ) + 2E(s p dp) 



(A.6) 



In our case, E(s 2 ) can be used to estimate E(s p 2 ). In the case 
where domain T> is an harmonic domain, p indexes harmonic 
modes ((, m), and E(s p 2 ) is a term of the CMB power spectrum. 
In our needlet approach, E(s p 2 ) is also directly connected to the 
CMB power spectrum. For this reason, it is important to char- 
acterise in the best way we can the "noise bias" E(d p 2 ) and the 
covariance of the error with the CMB, E(s p d p ). 

The ILC being constructed so that the response to the signal 
of interest is unity, only the filtered noise term contributes to the 
error d„, which can then be written as: 



dp = 



a'R: 



a' R x l a 

a' [R x + A,]- 1 

i — n P 

a' [R x + AJ- 1 a 



(A.7) 



(A.8) 



where A A is a small perturbation to R A . We use the first order 
expansion: 

[R + AJ" 1 * R" 1 - R" 1 A x R- 1 



which yields 

a' [R- 1 - RIJ'A.R; 1 ] n p 
dp ~ a' [R-i - R^R-i] a 
Writing: 

1 



(A.9) 



(A. 10) 



1 



1 



a' [R; 1 - FVA^R-i] a 



where e is 

a'R: l A r Rz l a 



a'R- x l a (1 
1 



a'R- 1 a 



e) 
(1 + e) 



a'Rz l a 



we get 



d p ^ 



a' [R; 1 - R^R; 1 ] n f 



(1 + 6) 



a'R- x l a 

Keeping only first order terms in (A x ) yields: 

_ a'R- x l np a'R- l k x R- l n p 
" ~ a'R-M 



(All) 



(A 12) 



a'R- l a 

[a'R x l np}[a'R x l A x R- l a] 

[a'R-'a] 2 



(A. 13) 



The first term on the right hand side, proportional to n, is the 
"ideal" ILC error, i.e. the error we would get if we knew per- 
fectly the "true" covariance matrix R A of the observations. The 
second and third terms, proportional to A x , are corrections due 
to the fact that this covariance matrix is actually obtained empir- 
ically from the observation themselves. 

From equations A. 3 and A.4 we can write A v in the form: 



A x = 5 s aa! + A„ + C 



(A 14) 
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where 



A„ - R„ - R„ 



(A.15) 
(A. 16) 

(A. 17) 



These three quantities correspond respectively, in pixel (or 
mode, or needlet coefficient) p, to the uncertainty in CMB vari- 
ance estimates due to "cosmic" (or sample) variance, to the er- 
ror in the estimation of the "noise" covariance matrix alone (if 
maps of noise+foregrounds alone were available), and to a cross 
term, originating from the empirical covariance between CMB 
and noise due the finite sample size N p 



The two last terms (small correction terms) in equation A. 13 



being proportional to A A , can be decomposed each into the sum 
of three terms, proportional to (S s aa'), A„, and C respectively. 
The signal and noise realisations enter the (6 s aa l ) term as prod- 
ucts of terms of the form (aa's q s q ti p ) only, the A„ term as prod- 
ucts of terms of the form (n q n' n p ). On the contrary, signal and 

noise realisations enter the C term as the product of terms in the 
form (as q ii' rip), i.e. products of s and second power of n. Index 
q runs over domain D. 

Assuming that s and n are centred variables, the mean value 
of the error is immediately seen to vanish: 

E(^) = (A. 18) 

The main contribution to the variance comes from the first term 



on the right hand side of eq. A. 1 3 The second and third terms are 
small corrections to this variance estimate, so that to first order, 
we get: 



E(^) - 



a'R-tRnRfa 

[a'R- l a] 2 



(A. 19) 



Recalling that R A = R„ + cr 2 aa', where a 2 is the variance of the 
CMB, and making use of the inversion formula: 



[R„ + (r 2 aa'l l = R„ 



we finally obtain: 



1 +o- 2 a'R- l a 



E(4) 



1 



[a'R-'a]- 



(A.20) 



(A.21) 



The most interesting terms are those connecting the error to 
the signal of interest, E(s p d p ), which is necessary to compute the 
power spectrum of the output map according to A.6| 8 



As mentioned previously, under the assumption that the sig- 
nal of interest is not correlated to the noise and the foregrounds, 
the first term (main term) of the r.h.s. of equation A. 13 does not 



give rise to multiplicative errors (or correlation of d p with s p ). 
Similarly, the corrective term proportional to 5 s aa' + A„, mul- 
tiplied by s, gives terms which contain an odd power of s and 
an odd power of n, and does not give rise to correlations. This 



8 We warn the reader that some authors fail to make a clear distinc- 
tion between the statistical (ensemble average) correlation, which is a 
deterministic quantity, and the "empirical correlations" computed, as- 
suming some kind of ergodicity, as averages over finite sets of samples 
as in equation A. 3 



assumption is correct, to excellent accuracy, when the signal of 
interest is CMB anisotropiesrlWe are left with: 



E(spd p ) = E 



y s p s q [« ,R ;S] [«'R;'(«,«' + an q )R- x l a\ 



z 

V q 



[a'R-^aY 
s q a'R^ (n q a' +an t q )R~ 1 n p y 



N n 



a'R: l a 



(A.22) 



Multiplying the numerator and denominator of the second term 
by a'R^'fl and expanding numerators, two terms cancel and two 
remain, f we assume in addition that signal and/or noise coeffi- 
cients are independent, i.e. E{s p s q ) = cr 2 5 qp and/or E(n p n') - 
Rn5 qp , only the pp term is non vanishing, and we get 



a ( 



E(s p d p ) 



N„ 



E ((a'R- 1 n p ) 2 ) - (a'R» E (^R; 1 ^)) 

[«'R>] 



(A.23) 



We compute 



E ((a'R; 1 »,)*) = a'R^RnR^a 

= a'R; 1 [R. v -o->a']R; 1 a 

= [«'R; 1 a][l-<r?a'R; 1 a] 
and 

E^R; 1 ^) = Tr^'R,,) 

= Tr^'fR,-^']) 
= Tr(ld)~cr 2 Tr(R;W) 
= m -cr 2 (a'R; l a) 



(A.24) 



(A.25) 



where m is the number of channels used for the ILC (here, 5 
WMAP channels + 1 IRIS map, for a total of 6). Substituting 
the results of equations |A.24 and A.25 into equation A.23 
get the simple final result: 



we 



E(s p d p ) = 



o-2(l -m) 



(A.26) 



The error in the reconstructed CMB map comprises a term pro- 
portional (on average) to the CMB. In our application, m - 6 



and N,, 



1024, so that if indeed all needlet coefficients were 



independent, the amplitude of the effect should be E(s p d p ) =* 
5 x 10 ~ 3 <t 2 , i.e. a bias of about half a percent in the CMB recon- 
struction. 

A.2. A geometric interpretation 

Although allowing the statistical derivation of the (anti-) correla- 
tion of the reconstruction error with the original CMB, the above 
derivation is not very illuminating about the mechanism giving 
rise to this CMB power loss. A geometrical reasoning gives bet- 
ter insight on what is actually going on. 

For a given data set, the ILC works on one single realisa- 
tion of all random fields. For an independent implementation of 



9 Certainly the CMB is not correlated to galactic components. Small 
correlations with large scale structure, and hence with SZ effect and 
emission from outer galaxies, may exist because of the integrated 
Sachs-Wolfe effect. We neglect this effect in the present discussion. 
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the ILC on N p pixels (or modes, or needlet coefficients) of the 
observations, each data set is represented by a vector in a N p - 
dimensional vector space W. The CMB s, the observation jc, for 
each channel i, and each of the noise realisations «, (including 
foregrounds) are elements of W. 

The collection of vectors «,- defines an m-dimensional sub- 
space V of W. This is true irrespective of the nature of the fore- 
grounds: indeed, although in principle vectors «, could be lin- 
early dependent, this happens in practice with vanishing proba- 
bility (in particular if the observations are noisy). 

Vector space W can thus be decomposed in two orthogonal 
subspaces, U and V, where V is the ra-dimensional sub-space 
spanned by all vectors «,, and U - V 1 " is a (N p - m) dimensional 
vector space. The CMB itself can be decomposed into two com- 
ponents, one in U, and one in V: 



s = Su + s v 



(A.27) 



where Sjj is the orthogonal projection of s onto U, and Sy its 
orthogonal projection onto V. 



U(dim 




V(dim.=m) 



Fig. A.l. Geometric illustration of the ILC bias. The true CMB 
(blue) and the output CMB (red) are elements of W - U © V. 
The difference between the two can be decomposed into the sum 
of two elements of V: a bias (green), which is an element of K, 
and a noise contribution (black), which is an element of H, the 
orthogonal of K in V. 

What happens when the ILC is made is the following: we 
look for weights w, for all channels, that minimise the variance 
of the reconstructed map, i.e. minimise the norm of vector T - 
2, WjXi, under the constraint that £, w, = 1. We have: 



s - Su + \Sy + 



Su + Sy 



I 



Willi 



(A.28) 



where the first term is a vector of U and the second term a vector 
of V, and the second line of the equation deftness^ and Ty. Since 
these two subspaces of W are orthogonal, the norm ofT is the 
sum of the norms of the two vectors Ty and Ty. The norm ofT 
thus depends on w, only through the norm of the projection ofT 
on subspace V. 

The noise contribution to T appears as a linear combination 
of vectors «,. For varying values of w t such that 2,W; = 1, 
this linear combination spans an affine subspace S of V. S is 
of dimension m - 1 (an hyperplane). Defining K as the vector 
subspace of V spanned by linear combinations £,- w,jc, such that 
Yji Wi - 0, we obtain S as: 



S =p + K 

where p is any element of S . 



(A.29) 



We note that the vector subspace K depends only on noise 
realisations, and not on s nor on the final ILC weights (the lat- 
ter only defining a single point on S - and on K by orthogonal 
projection). Hence, the direction of the one-dimensional vector 
subspace H of V orthogonal to K is also independent of s and of 
the final ILC weights. 

For any element Tv = Sy + 2j w, rij of affine space S the norm 
of Ty is the sum of the norms of its projections onto K and H. 
Allowing weights w, to vary, vector Ty spans S , and hence only 
the norm of the projection of Ty onto K varies (and not its pro- 
jection on H). The minimum is reached when the projection of 
of ~sy onto K vanishes. When this happens, the ILC has cancelled 
completely the linear combination of projections of s and «, onto 
the m - 1 -dimensional space K, and left untouched the projec- 
tions of s and n, onto the N p - m + 1 dimensional vector space 
U@H. 

Assuming the CMB to be Gaussian and statistically 
isotropic, its coefficients in any orthogonal basis are Gaussian 
distributed random variables with variance cr 2 /N p (since the sum 
must have total variance o 2 ). It follows straightforwardly that the 
correlation of the recovered CMB map with the input CMB map 
is (N p - m + 1)/N P , and that the "bias" is due to the loss of m - 1 
modes of the original CMB, which have been unlucky enough to 
"live" in the (m - 1) dimensional space K. 

A.3. Comment on coefficient independence 



The above derivation in Section |ATTj assumes the independence 
of coefficients n p and/or of coefficients s p , i.e. E{n p n') = R„5 ?p 
and/or EispSq) = o- 2 s 6 qp . When this assumption does not hold, 
we have: 

-^Z4v ? Kfl;\)) (A.30) 

Assuming that the noise and the CMB are independent, we have: 

I E((q'FL> p )(q'FL:S)) 
E (Vp) = Wp h E(v^ F^^ 

-^Z E W E K R ^) (A - 31) 

where 

E((« r R; 1 H p )(a'R; 1 «,)) = a'R- l E(n p ri q ) R~ x a (A.32) 

and 

E^R- 1 *,) = Tr(R;'E (n p n' q )) (A.33) 

When p and q index needlet coefficients as in the present 
work, we have: 



2^+1 



E { s p s i) = X ~n — h2 e c e p c(cos6 qp ) 



(A.34) 



where s p and s q are needlet coefficients of the CMB map, eval- 
uated at two different points p and q, Q is the angular power 
spectrum of the CMB, 9 qp is the angle between q and p, and 
N tot is the total number of pixels of the needlet coefficient map. 
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For noise maps (including foregrounds), which are not station- 
ary Gaussian random fields on the sphere, the analogous formula 
is just an approximation, which can be written as: 

E (n p n' q ) - Y 2 -^h] R n {€) P ( (cos qp ) (A.35) 

where n p and n q are needlet coefficients of all noise maps, eval- 
uated at two different points p and q, R„({) is the covariance of 
the noise needlet coefficients (an rax ra matrix for each €), and 
qp and N tot are defined as above. 

Assuming that neither the noise level, nor the CMB power, 
do vary much over the spectral window hi, R„(£) and Q are 
approximately independent of (, and can be taken out of the in- 
tegral. We get: 

E( SpS q)=C e k(e qp ) (A.36) 

and 

E(n p n' q ) = R n k(6 qp ) (A.37) 

with 

**M = X ^r^A/^(cos 9 qp ) (A.38) 



Hence, plugging this result, together with equations A. 32 and 
A. 33 into equation |A.31| we get: 



(1-ra) 



E( Sp d p ) = K —— -!- C[ Y k 2 (8 qp ) (A.39) 

"p q 

Finally, noting that <x 2 - Cik(0), we get 

EM,) = %^ j^j 2 ^ ( v ( A -40) 



Equation A. 40 is the exact same as equation A. 26 except for 
a coefficient, which measures the correlation between signal and 
noise coefficients p and coefficients q in domain D. In particular, 
the result is, again, independent of R„. 

Hence, we define an effective number of modes, N p = N p /a, 
where 

a = ~1W~ (A ' 41) 

and we get 

0-1(1 -m) 
E(d pSp ) = p (A.42) 

Nf 

An approximation (and upper bound) for a is easily obtained in 
the special case where he is a square spectral window, and when 
domains T> over which the ILC is implemented are small regions 
of the sky, so that h\ - 1, and P^(cos 6 qp ) =* 1. In this case, we 
have: 



^(M-T?-((^ax + l) 2 -4n) 



(A.43) 



N tot 
and we have 

<-|^((W + l) 2 -4») (A-44) 

We note that ((^ max + l) 2 - l^in) simply is the number of har- 
monic coefficients selected by the spectral window hg, and that 
N p /Ni t = / s ky is a coefficient which takes into account the effect 
of partial sky coverage for the (local) calculation of the statis- 
tics of the data set, well in line with N'f being understood as an 
effective number of modes. 



